### Read in the libraries
library(tidyverse)
library(rstan)
library(modelsummary)

### Read in the data
data <- read.csv("post96_data_trimmed.csv")

# Remove cases with zero variance 
cases_with_variance <- post96 %>%
  group_by(c.casenr) %>%
  summarize(thevar = var(v.vote, na.rm = TRUE)) %>%
  filter(thevar > 0) %>%
  pull(c.casenr)

post96 <- post96 %>%
  filter(c.casenr %in% cases_with_variance)

# Make justice data for later analysis
post96_judges <- post96 %>%
  distinct(j.name.y, .keep_all=T) %>% 
  dplyr::select(j.judgenr, 
         j.name.y,
         j.appointinggov, 
         j.independentappointment,
         j.career, 
         j.attorney, 
         j.governmentattorney, 
         j.legalacademic, 
         j.legislationdep, 
         j.publicprosecutor,
         j.interim_only,
         j.sex,
         j.earlierjudge, 
         j.independentappointment,
         j.bornoslo) %>% 
  mutate(judge = j.name.y, 
         j.governmentadministration = ifelse(
         j.legislationdep == 1|j.governmentattorney == 1|j.publicprosecutor == 1, 1, 0),
         j.career_government = ifelse(j.career == "Government adm.", 1, 0),
         j.career_attorney = ifelse(j.career == "Attorney", 1, 0),
         j.career_judge = ifelse(j.career == "Judge", 1, 0),
         j.career_academic = ifelse(j.career == "Academic", 1, 0))

# For descriptive table 
judges <- post96_judges %>% 
  dplyr::select(j.appointinggov, 
         j.independentappointment,
         j.career_government, j.career_attorney,
         j.career_judge, j.career_academic, 
         j.governmentadministration, j.attorney,
         j.legalacademic, j.earlierjudge, 
         j.legislationdep, j.governmentattorney, j.publicprosecutor,
         j.sex, j.bornoslo, j.interim_only)


# For Stan
votes <- post96$v.vote
N <- length(votes)
j <- as.numeric(factor(as.character(post96$j.name.y)))
jlevels <- levels(factor(as.character(post96$j.name.y)))
k <- as.numeric(factor(as.character(post96$c.casenr)))
klevels <- levels(factor(as.character(post96$c.casenr)))
J <- max(j)
K <- max(k)

leftanchor <- which(jlevels == "Ingse Stabel")  
rightanchor <- which(jlevels == "Jens Edvin A. Skoghøy")

# Now do a more structured model. To achieve these, we need to
# collapse the case-level parameters to a data frame
casedf <- post96 %>%
    dplyr::select(c.casenr,
                  civil = c.civil.x,
                  echr = c.echr,
                  eea = c.eea_eu,
                  economy = c.economy,
                  constitutional = c.constitutional.reference,
                  enlarged = c.grandplenary,
                  criminal = c.criminal,
                  sentencing = c.criminalsentencing.kw,
                  pro_public_civil = pro_government_civil,
                  pro_public_echr = pro_government_echr,
                  pro_public_eu = pro_government_eu,
                  pro_public_economy = pro_government_economy,
                  pro_public_constitution = pro_government_constitutional,
                  pro_public_enlarged = pro_government_enlarged_panels,
                  pro_public_criminal = pro_government_criminal,
                  pro_prosecutor_sentencing = pro_government_sentencing) %>%
    group_by(c.casenr) %>%
    summarize_all(mean) %>%
    replace(is.na(.), 0)

## Ensure the correct ordering
casedf$c.casenr <- factor(casedf$c.casenr,
                          levels = klevels,
                          ordered = TRUE)
casedf <- casedf %>%
    arrange(c.casenr)

XX <- model.matrix(~1 + pro_public_civil  
                  + pro_public_echr + pro_public_eu + pro_public_constitution
                  + pro_public_economy + pro_public_enlarged 
                  + pro_public_criminal + pro_prosecutor_sentencing,
                  data = casedf)

nvarsXX <- ncol(XX)

hier_data_large <- list(N = N,
                  K = K,
                  J = J,
                  X = XX,
                  nvars = nvarsXX,
                  j = j,
                  k = k,
                  y = votes,
                  leftanchor = leftanchor,
                  rightanchor = rightanchor)



## ----stanmodel-----------------------------------------------------------

stanirthier <- '

data {
  int<lower=1> J; //Legislators
  int<lower=1> K; //Proposals/bills
  int<lower=1> N; //no. of observations
  int<lower=1, upper=J> j[N]; //Legislator for observation n
  int<lower=1, upper=K> k[N]; //proposal for observation n
  int<lower=0, upper=1> y[N]; //vote of observation n
  int<lower=1, upper=J> leftanchor; // anchored to -1
  int<lower=1, upper=J> rightanchor; // anchored to +1
  int<lower=1> nvars; // number of explanatory coefficients
  matrix[K, nvars] X; // model matrix
}
parameters {
  vector[K] alpha;
  vector[K] beta;
  vector[J] theta;
  vector[J] mu_theta;
  real<lower=0> sigma; // standard deviation on case regression
  vector [nvars] coefs;
}
model {
  alpha ~ normal(0,10); 
  theta ~ normal(0, 1);
  theta[rightanchor] ~  normal(1, .01);  //constraints
  theta[leftanchor] ~ normal(-1, .01);

  coefs[1] ~ cauchy(0, 10);
  for (i in 2:nvars) {
    coefs[i] ~ cauchy(0, 2.5);
  }
  beta ~ normal(X * coefs, sigma);
  
  y ~ bernoulli_logit(rows_dot_product(theta[j], beta[k]) - alpha[k]);    
}
generated quantities {
  real<lower=0, upper=1> pcp;
  vector[N] cp;
  int<lower=0, upper=1> y_new[N];

  for (n in 1:N) {
      y_new[n] = bernoulli_logit_rng(theta[j[n]] * beta[k[n]] - alpha[k[n]]);
      cp[n] = y_new[n] == y[n]; 
  }
  pcp = mean(cp[]);
}

'



## ----loadprecomputed-----------------------------------------------------
if (file.exists("./stan_fits.RData")) {
    load("./stan_fits.RData")
} else {
    hier.fit.large <- stan(model_code = stanirthier,
                           pars = "cp",
                           include = FALSE, # exclude cp
                           data = hier_data_large,
                           iter = iter_incl_warmup,
                           chains = 2,
                           thin = 5,
                           init = "random",
                           refresh = 10,
                           cores = 2, seed = 1234,
                           control = list(max_treedepth = 15))
    
    save(hier.fit.large,
         file = "stan_fits.RData")
}




## ----getfit--------------------------------------------------------------
## Model fit statistics

## Calculate the GMP
calcGMP <- function(ystar, y) {
    pijt_1<-log(ystar)
    pijt_0<-log(1-ystar)

    pijt_1<-pijt_1*(y==1)
    pijt_0<-pijt_0*(y==0)

    my.logLik<-sum(pijt_0,pijt_1,na.rm=T)

    my.A<-length(which(!is.na(y)))

    return(exp(my.logLik/my.A))
}

predprobs <- summary(hier.fit.large, pars = "y_new")$summary[,"mean"]
gmp.hier.large<-calcGMP(predprobs, hier_data_large$y)
pcp.large<-summary(hier.fit.large, pars = "pcp")$summary[,"mean"]

predprobs0 <- mean(votes)
gmp0<-calcGMP(predprobs0, hier_data_large$y)

PCP<- c(predprobs0, pcp.large)
GMP<- c(gmp0, gmp.hier.large)
model.fit<-as.data.frame(cbind(PCP, GMP)) %>% 
  mutate(Model = as.character(c("Null", "Hierarchical IRT"))) %>% 
  dplyr::select(Model, PCP, GMP)



## ----judgeplots----------------------------------------------------------

hier_judge<- summary(hier.fit.large, pars = "theta")$summary %>%
  as.data.frame() %>%
  dplyr::select(mean = mean,
         median = `50%`,
         lo = `2.5%`,
         hi = `97.5%`,
         sd = `sd`,
         rhat = `Rhat`) %>%
  mutate(judge = 1:n())


hier_judge$judge <- jlevels

# merge ideal points with justice attributes (also text vars for ideal point graphs)

ideal_judge<-hier_judge %>% 
  inner_join(post96_judges, by = "judge") %>%
  mutate(j.appointinggov.text = ifelse(j.appointinggov == 1,
                                       "Social democratic",
                                       "Conservative"),
         j.independentappointment.text = ifelse(j.independentappointment == 1,
                                                "Post-reform",
                                                "Pre-reform")) 

ggplot(ideal_judge, aes(x = reorder(judge, (median)), y = median,
                         ymin = lo, ymax = hi)) +
  geom_pointrange(aes(shape = factor(j.appointinggov.text), color = factor(j.independentappointment.text)), 
                  show.legend = TRUE) +
  scale_shape_manual(values = c(0, 17, 2, 1)) +
  scale_color_manual(values=c('black','dark grey')) +
  coord_flip() +
  theme_minimal() +
  labs(y="Ideal Point", x= "") +
  theme(text = element_text(size=10)) +
  theme(legend.position = "bottom",
        legend.title = element_blank())



## ----appointments--------------------------------------------------------

lm_appointment <- lm(median ~ j.appointinggov + j.independentappointment, data = ideal_judge) 

lm_appointment2 <- lm(median ~ j.appointinggov * j.independentappointment, data = ideal_judge)

lm_career1<-lm(median ~ j.career,
               data = ideal_judge)

### Export regression tables
cm <- c('j.appointinggov' = 'Appointed under Soc-Dem govt',
        'j.independentappointment' = 'Appointed after 2002 reform',
        'j.appointinggov:j.independentappointment' = 'Soc.-Dem appointee post 2002 reform',
        'j.careerAttorney' = 'Formerly lawyer in private practice',
        'j.careerGovernment adm.' = 'Former civil servant',
        'j.careerJudge' = 'Career judge')

cat(reg_tab <- msummary(list(lm_appointment, lm_appointment2,
              lm_career1),
          coef_map = cm,
          stars = c('*' = .05, '**' = .01)) %>%
    knit_latex())

### Or alternately just print(reg_tab)


